Scale Invariance in disordered systems: the example of the Random Field Ising Model 
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We show by numerical simulations that the correlation function of the random field Ising model 
(RFIM) in the critical region in three dimensions has very strong fluctuations and that in a finite 
volume the correlation length is not self-averaging. This is due to the formation of a bound state 
in the underlying field theory. We argue that this non perturbative phenomenon is not particular 
to the RFIM in 3-d. It is generic for disordered systems in two dimensions and may also happen in 
other three dimensional disordered systems. 
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Ferromagnetic systems in a random fields (RFIM) have 
a puzzling behaviour. It is already known that perturba- 
tive renormalization group (PRG) leads to the prediction 
of dimensional reduction that is known to give wrong re- 
sults for the RFIM, and this although RFIM is, together 
with branched polymers 0, the only case where PRG 
can be analyzed to all orders of perturbation theory jl],U . 
Several attempts were made to explain this discrepancy 



One suggested possibility is the appearance of a non 
perturbative phenomenon (i.e. the formation of a bound 
state in the underlying field theory) which invalidates the 
perturbative analysis. There are several consequences of 
this. It is not any more obvious that naive scaling argu- 
ments apply in the scaling region as it is usually assumed, 
or if and how they should be modified. We expect also 
these non perturbative phenomena, related to the forma- 
tion of bound states, to invalidate the predictions of per- 
turbative renormalization group for the critical behaviour 
of disordered systems (calculation of critical exponents, 
proof of universality classes etc.) A new examination of 
these questions is required. 

We claim that the appearance of bound states is not 
peculiar to the 3-dimensional RFIM (what is peculiar to 
the RFIM is the possibility to perform the very large 
scale simulations needed in order to observe these non 
perturbative phenomena). Kardar et al. [[)) have already 
noticed that averaging over the disorder generates at- 
tractive interactions among replicas that produce bound 
states in several two dimensional systems. In the for- 
mation of bound states there is competition between the 
strength of the attractive forces and the size of the avail- 
able phase space. In two dimensions phase space is small 
and even a very small attraction wins, so bound states are 
generic. In three dimensions formation of bound states 
depends on the strength of the attraction. In the case 
of the RFIM Brezin and De Dominicis || have already 
found that an instability appears in the kernel of the 
Bethe-Salpeter equation, leading to the formation of a 



bound state in less than 4 dimensions. Our results con- 
firm their analysis. 

The relevance of the attractive interaction among repli- 
cas can also be seen in a case where the interaction be- 
tween replicas is repulsive. This is the case of diluted 
branched polymers (or lattice animals). There are not 
bound states in that case and perturbative renormal- 
ization group, i.e. dimensional reduction, has to work. 
Indeed dimensional reduction has recently been proven 
rigorously in this case |2(]] . 

In this paper we address the question of the existence of 
bound states and of the validity of the usual scaling laws 
by studying the 3-d RFIM by extensive numerical simu- 
lations at zero temperature. As many previous authors 
pr|- p^| , pT[ we have taken advantage of its equivalence at 
zero temperature with the maximum flow problem in a 
graph p0|, for which a very fast polynomial algorithm is 
known Jl8[ | . The simulations mentioned above provided 
evidence of an infinite correlation length at the phase 
transition through a finite size scaling analysis, but no di- 
rect measurement of a correlation length was performed. 
In the present paper we perform direct measurements 
of correlation lengths. We will see that the correlation 
functions have very strong fluctuations. We will show by 
a scaling analysis that these fluctuations are so impor- 
tant, that in any finite volume the correlation length in 
the critical region is not self-averaging (i.e. it is strongly 
sample dependent). 

Let us remind that the Hamiltonian of the RFIM is of 
the form H = —JJ2<ij> a i a j ~ J2i hi^i- As usually 
J2<ij> runs over neighbouring sites of a cubic lattice 
with coordinates x,y, z — 1, • • • , L and hi are indepen- 

2 

dent random Gaussian variables with variance h t = 1. 
It is well known that there is a zero temperature phase 
transition in this model. For J > J c the spins are fer- 
romagnetically ordered, while for J < J c we get a dis- 
ordered phase where a large number of spins are locally 
aligned with the random external field. 
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FIG. 1. Masses m(k) 1/v (with v = 1.5 ) as a function of 
J, for L = 160. k = 1 is the lowest curve, while k — 3 is the 
highest. 

In order to measure correlation length at zero temper- 
ature we proceed as follows. We consider cubic lattices 
with linear size L. For every sample of the random field, 
we considered two copies of the system with different 
boundary conditions. For the first copy a all the spins 
at the plane x = 1 are set to plus one. For the other 
copy t they are set to minus one. We let free boundary 
conditions for both copies at the other end of the lattice, 
i.e. x = L. We impose periodic boundary conditions on 
the two perpendicular directions. For each value of the 
ferromagnetic coupling J we find the two ground states 
<Ji and Ti . For every sample s of the random field we mea- 
sure the distance d s (x) = s(l — Qs( x )) where the overlap 
is given by 
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which depends on the sample s; (we remind that 
because of the boundary conditions at i = 1, 
cr(l, 2/0, zo)t(1, Vo, z o) = — 1 ) d s (x) measures the propor- 
tion of spins which have opposite values in the two ground 
states at distance x from the boundary. By construction 
d s (l) = 1. If J << J c d s (x) will rapidly decrease with 
x, while for J » J c it will asymptotically remain con- 
stant. It is known from lattice gauge theories |2l| that 
one efficient way to measure correlation lengths, is to im- 
pose boundary conditions on an observable and measure 
the variation of the observable with the distance from the 
boundary. From the behaviour of d s (x) measured numer- 
ically we are going to measure the correlation length and 
draw conclusions on the underlying field theory. 

We simulated 1200 samples for L = 160, 2000 samples 
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FIG. 2. Masse-ratios as a function of J, for L = 160. Cir- 
cles represent m(2)/m(l) and squares represent m(3)/m(l) 

for L — 120 and 3500 samples for L — 80 for different 
values of the ferromagnetic coupling J < J c . From the 
previous simulations p^- p^Jl7) | it is known that J c ~ .44. 
We found very strong sample to sample fluctuations of 
d s (x), even for systems of sizes 160 3 . In order to extrapo- 
late to the infinite volume limit, we studied the behaviour 
of d^ k '(x) — d s {x) k for k — 1,2,3, (as usually the bar de- 
notes average over the random fields), i.e. the average of 
the /c'th power of d s (x) over the samples. 

For J < J c we expect an exponential decay in x of 
d^{x) ~ exp (— m(k)x). The masses m(k) are inverse 
correlation lengths. We first fitted our data with the 
formula 



d(k) = exp(-m(fc)x) 1 



(2) 



We expect m(l) to vanish at the critical value J = J c . 
According to finite size scaling m(l) ~ ( J c — J)"f((J c — 
J)L 1 / V ). We found that the scaling function / is a con- 
stant with a good approximation, i.e. a very small L 
dependance of m(l), and 1.3 < v < 1.9, compatible with 
previous results |l4|,[l^,[llj . We have obtained a better 
determination of v in another set of simulations, (with 
different boundary conditions) that we will publish else- 
where. For a given value of J c the statistical errors in v 
are very small. The large uncertainty is due to subdomi- 
nant corrections and to the extreme dependance of v on 
the value of J c . For J « J c (in practice for J < .43) 
we found a(l) = in agreement with dimensional argu- 
ments, while for J ~ J c a(l) = .27 ± .05. 

In figure 1 we plot m{k) 1 ' v (with v = 1.5 ) as a func- 
tion of J for L = 160. Not too far from J c , to(I) 1 / 1 ' is 
linear in J to a very good approximation, to(2) 1 /" and 
m(3) 1//l/ have a more complex behaviour. What is very 
surprizing is that for J > .425, m(2) and m(3) collapse 
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FIG. 3. Probability distribution (histogramme) of d(x) 
for J = .436, L = 160 and two values of x, x = 55 (dashed 
line) and x = 110 (continuous line). 

with m(l). This is best illustrated in figure 2, where 
we plot the ratios m(2)/m(l) and m(3)/m(l) as func- 
tion of J for L = 160. Almost identical results were 
obtained for L = 120 and L — 80. For small enough 
ferromagnetic coupling m(k)/m(l) — k, i.e. the masses 
are self-averaging, as expected in perturbative field the- 
ory |nj], while close to the critical point m(fc)/m(l) = 1. 
It is obvious that < d s (x) < 1 from the very definition 
of d s (x). It follows that d^(x) < d^(x) and therefore 
m(k) > to(1); m(k) = m(l) amounts to the maximum 
possible violation of self-averaging of d s (x). This is also 
illustrated in figure 3, where we plot the probability dis- 
tribution (histogramme) H(d) of d s (x) over the random 
field samples, for L — 160, J = .436 and and two val- 
ues of x, x = 55 (dashed line) and x — 110 (continu- 
ous line). Theses histogrammes separate into two peaks, 
one around d = and the other around d = 1. As we 
vary x the shape of the histogramme remains the same, 
only the relative height of the two peaks changes. For 
smaller but not too small J's, in the region of J where 
still m(k) ~ to(1), H(d) has also contributions outside 
the two peaks. The peak around d = 1 is necessary 
for m(k) — to(1), otherwise d^ k '(x) would decay with x 
faster than d^'(x). 

For any finite x there is also the possibility of a su- 
perposition of exponentials in d^ (x) . We have therefore 
tried another fit to d^(x), 

d (2 \x) = Cl d {1 \x) +c 2 {d {1) {x)) 2 +c 3 exp(-m'a;) (3) 

where m is a new masse. It turned out, particularly 
for intermediate values of J, that this is a much better 
fit. We found that no superposition of exponentials is 
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FIG. 4. Scaling plot of cl versus mL (see text). Circles 
correspond to L = 160, squares to L = 120, and stars to 
L — 80. 

needed in the fit of d^(x). We found that when C3 is 
7^ 0, m > 2m(l), i.e. this term contributes only for 
small x. 

g?W (x) has the slowest decay in x in equation a. There- 
fore ci measures the violation of self-averaging of d at 
large x. For small J, c\ ~ C3 ~ i.e. the masse is self- 
averaging, while for J ~ J c , c% ~ C3 ~ 0, c\ ~ 1. In order 
to study the crossover between the two regimes, we tried 
the finite size scaling ansatz 

ci - g{{J c -J)l}l v ) = g{Lm) (4) 

where m = m(l). In the previous equation we used the 
finite size scaling of m, i.e. mL — f((J c — J)L X ^ V ) to 
change variables from (J c — J)!?-! 1 ' to mL. As is seen 
in figure 3 where c\ is plotted versus mL this ansatz 
works rather well. In the scaling region, i.e. for small 
to, ci ~ 1 and d^(x) ~ exp — mx, i.e. m(2) = m(l) 
and we observe the maximum possible violation of self- 
averaging, as we have seen before. Outside the scaling 
region (i.e. for to > .075) we found that for constant J, 
ci ~ 1/L 2 (we remind that in this region m depends on 
J and not on L ) and d^ 2 \x) ~ exp— 2mx. This means 
that the violation of self-averaging of d vanishes as 1/L 2 . 
From the very definition of d s (x), if d s ^ 0, d s > 1/L 2 . It 
follows that d^(x) > d(x)/L 2 7 i.e. c\ ~ 1/L 2 amounts to 
the smallest possible violation of self-averaging. It is, we 
think, very remarkable that in the scaling region we cross 
from the smallest possible violation of self-averaging, i.e. 
ci ~ 1/L 2 , over to the maximum possible violation of 
self-averaging, i.e. c\ ~ 1. 

The quantities d'" have a simple correspondence in 
field theory at finite temperature T . Let us consider an 
infinite system at finite temperature. We study the field 
theoretic equivalent to the quantity q s (x) defined in equa- 
tion p\ and its correlations. To the Ising spins a and r 
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correspond the fields <p(x) and tp(x). Their vacuum ex- 
pectation (i.e. their space average) has been substracted 
from <p(x) and ip(x) so that their correlation, for exam- 
ple (4>(x)(j)(y)) decay with their distance \S — y\. The 
connected correlation function is defined as 

G s (r,f ') - (0(r)^(r)0(r >(r ')> (5) 

The subscript s represents the dependence on the ran- 
dom field sample. We can also define G^ k \f — r ') the 
average over the samples of the fcth moment of G s (r,f') 
(averaging over the random fields restores translation in- 
variance), G^ k \r—r') = (G s (r,r')) k . These correlations 
G^ k \f — r ') can be simply defined if, as usually, we in- 
troduce n replicas (at the end n — ► 0) to average over the 
random fields. In the high temperature phase and for 
large \r-f'\, G^ k \r-r') ~ cxp 'D If W e sum 

over the "transverse" y and z components of the vectors 
r and f ', G' {k \r x - r' x ) = J r r> f , r , G^ k \r- f '), then 

G'W(r x -0 ~ g^M^^ i.e. G^W^-O and 

G^ k \r — r ') have the same exponential behaviour, only 
the power prefactor is different. If we define 

y,z 

where <j> and tp are the fields corresponding to two iden- 
tical copies of the system, and the "composite" oper- 
ators A ala2 - a " (x) = Q al {x)Q a2 (x)---Q ak (x), where 
a 1 ,a 2 ■ ■ ■ a k are k replica indices all different from each 
other, 

G {k \x-y) = {A ala2 - a \x)A ala2 - a \y)) (7) 
We can insert a complete set of states r in equation [?] 
G {k) {x-y)=^2{A ala2 - a \x)\r){r\A ala2 - ak {y)) (8) 

r 

Let's call ro(fc) the lowest mass state giving a non zero 
contribution to the previous sum and /x(fc) its masse. For 
large |x — y\, G ^(x — y) ~ exp — (fi(k)x). In perturba- 
tive field theory these lowest masse states ro(k) are those 
created by the 2k fields 4> a , tp al , Qt , i- e -j 

/i(fc) = 2kn, where /i is the masse of the field 4>. If we 
find fi(k) ^ 2/c/i perturbation theory brakes down. 

If universality is valid along the transition line, down 
to zero temperature, as it is usually assumed, the previ- 
ous arguments are also valid at zero temperature. First 
remark that for J < J c , and for x » 1, q s {x) = 1, so 
the "connected" part of q s (x) is q s . c (x) — q s {x) — 1 = 
—2d s (x). It follows from the previous that the long dis- 
tance behaviour of and G ^ (x — y) are the same. We 
found in our simulations that, near the critical region, 
m{k) ^ fcm(l), i.e. there exists a new state tq, coupled 
to all the composite operators A ^" 2 '" ^ (we have only 



measured the correlation functions for k = 1,2,3). This 
means that there exists a bound state and this explains 
why perturbation theory breaks down in the RFIM. 

A preliminary version of this work was presented by 
one of the authors (NS) at the Helsinki Workshop Disor- 
dered systems at low temperatures and their topological 
properties (January 2002). 
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